Forces and atomic relaxations in the pSIC approach with ultrasoft pseudopotentials. 
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We present the scheme that allows for efficient calculations of forces in the framework of pseudopo- 
tential self-interaction corrected (pSIC) formulation of the density functional theory. The scheme 
works with norm conserving and also with ultrasoft pseudopotentials and has been implemented in 
the plane-wave basis code quantum ESPRESSO. We have performed tests of the internal consistency 
of the derived expressions for forces considering ZnO and Ce02 crystals. Further, we have performed 
calculations of equilibrium geometry for LaTiOa, YTiOs, and LaMnOa perovskites and also for Re 
and Mn pairs in silicon. Comparison with standard DFT and DFT-I-U approaches shows that in the 
cases where spurious self-interaction matters, the pSIC approach predicts different geometry, very 
often closer to the experimental data. 

PACS numbers: 31.15.es, 61.50.Ah, 61.72.Bb, 61.72.S-, 71.15.-m 



I. INTRODUCTION 

Predictive power of the density functional theory^i 
mostly in its local density (LDA) and gradient corrected 
(GGA) flavors, is the main factor that has established 
this method as the standard approach in the materi- 
als science. For many electronic systems, it has be- 
come possible to predict very accurately the equilib- 
rium geometry, equation of state, relevant energetics, 
and further whole plethora of properties with astonish- 
ingly good accuracy. Unfortunately, all these approx- 
imations are plagued by the fact that functionals con- 
tain spurious self-interaction and the electronic states 
are typically too extended. Therefore, the reliable pre- 
dictions for systems with very localized electronic den- 
sity, so called strongly-correlated systems, require a com- 
putational scheme that cures the self-interaction prob- 
lem. In some approximate way, the self-interaction is 
partially removed in the DFT-I-U scheme^ which cor- 
rects the Coulomb potential within the localized states, 
such as d- and /-shells of atoms. There were also devel- 
oped methods with the exact exchange^'"— and the self- 
interaction correction (SIC)i^"— Perhaps the most simple 
among the DFT-I-SIC approaches, is the pseudopotential 
SIC (pSIC) scheme proposed by Filippetti et alr^ Its 
usefulness to reliably predict energetics has been widely 
proved in a variety of systems, to mention just a few such 
as transition metal oxides, manganites and cuprateSf^^ 
diluted magnetic semiconductors (DMSs) fi2>i^ strongly- 
correlated superconductors,— moleculesi^ and molecular 
junctions,''''' and many other as described in excellent re- 
view paper ^12 

Interestingly, the strongly-correlated systems exhibit 
very often strong deformations of the crystal lattice 
structures. The interesting and important examples in- 
clude Jahn- Teller distortions, relaxations around defects, 
atomic reconstructions at interfaces, lattice distortions 
due to magnetic interactions, surface reconstructions and 
local adjustment of atomic positions at surfaces due to 
the adsoption of atoms and molecules, and finally clusters 



of atoms in nanoparticles. It is obvious that the possibil- 
ity to calculate forces and stress tensor, in addition to the 
energy spectrum, consistently within the self-interaction 
free DFT scheme is very desirable. 

However, unfortunately, the full equations for forces in 
the pSIC method have not been set up yet and only an 
attempt to calculate forces, albeit in a very approximate 
form, has been performed in the paper by Filippetti and 
FiorentiniJ^ Even these simplified equations for forces 
have not been tested so far in any system. Only recently, 
a new variational pSIC approach,^* different than the 
original pSIC approach of Filippetti and Spaldin,'^ has 
been proposed. 

In this work, we provide a computational scheme that 
is based on the non- variational pSIC method,'^ imple- 
menting it into widely used quantum espresso codoi^ 
using the plane-wave basis and employing ultrasoft pseu- 
dopotentials (USPPs)j22 For this scheme, we also derive 
and implement the formulae for forces. It turns out that 
the procedure to calculate forces is similar to the one 
employed in the DFT+U method. 

The developed formalism is tested in a series of calcu- 
lations for various systems. We calculate internal strain 
parameter u for the wurzite ZnO and compare to the 
DFT-I-U results for the norm-conserving (NCPP) and the 
ultrasoft pseudopotentials. We perform tests also for the 
rare earth compound Ce02 with f valence shells. The 
relaxations of atomic positions in a cell are also tested 
for three chosen perovskites in distorted Pnma structure, 
namely LaTiOa, YTiOs, and LaMnOa. As a third test, 
we consider pairs of Mn and Re impurities in the silicon 
lattice, just addressing the problem of transition-metal 
ions pairing, that is so important for a relevant class of 
materials, namely the diluted magnetic semiconductors. 

The paper is organized as follows: the details of the 
implementation of the pSIC are given in section II, the 
full equations for forces are presented in section III, the 
illustrating implementations of the developed formalism 
are discussed in section IV, finally, the paper is summa- 
rized in section V. 
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II. IMPLEMENTATION OF THE PSIC 
METHOD FOR PLANE- WAVE BASIS 
COMPUTATIONAL SCHEME 

In this section we describe briefly all details necessary 
to implement the pSIC scheme, just to introduce unique 
notation necessary for section III. We follow closely for- 
mulation from the work by Filippetti et a/.,^^'^^ and col- 
lect here the most important equations. Note that the 
second paper— of the authors on this topic differs in some 
points from the first one,>ii mostly by setting additional 
simplifications which essentially do not affect accuracy 
but lead to a speed up of calculations. 

In the pSIC method, the Kohn-Sham equation for spin 
a orbitals (it implies the usage of a spin-polarized DFT 
approach) is corrected by the SIC potential Vsic 

[-V2 + Vpp + VSxc - Vsic] ICk) = <k ICk), (1) 

which is cast in the Kleinman-Bylander form,^^ and con- 
tains contributions from the all relevant local orbital po- 
tentials related to the local pseudo-orbitals 0i(r) (with 
index i describing lumped together angular momentum 
quantum numbers and position of the atom in the lattice) 
as follows 

^sic - / . — ^ — ■ 

i * 

The projection operators and the normalization in- 
tegrands Cf are defined 

7r(r) = ^^xcK(r)] U^), 

The V^jjf c potential is a sum of the Hartree potential 
and the exchange-correlation potential in a form that re- 
sults from the DFT functional used in the calculations. 
The V§xc potential is a functional of the local density 
nf (r) that is defined through the atomic pseudo-orbitals 

(r) and the occupation numbers pf 

nnr)^p7\Mrr, (3) 

^T.fnl.{rnM{^^\rni.)■ (4) 
nk 

The occupation numbers pf are obtained like in the 
DFT-fU scheme from the projection of the Kohn-Sham 
states ipnk onto the local atomic orbitals (f)i, and /^j^ are 
the Fermi-Dirac occupations. 

Note that if the pseudo-orbital functions do not de- 
pend on spin (as in a spin independent PF scheme used 
throughout this paper), the spin dependence of nf (r) en- 
ters only via the occupation numbers pf . 



It is important to perform orthonormalization of the 
local pseudo-orbital functions before using them in the 
above definition of pf, since it may change considerably 
the relations between the occupations of different atomic 
shells. This orthonormalization is not mandatory in the 
DFT+U method, since this scheme usually involves only 
one shell of given atom, d- or /-shell, but not the both. 

Further, the pSIC potential is scaled by one half for 
the relaxation contribution in the extended systems^i 

VSxcK] ^ I VSxcK]- (5) 

In general, the scaling coefficient is applied in this place 
to unify the bulk and molecular systems. ^^'^^ 

In order to simplify calculations, two approximations 
are made for the pSIC potential: 

1) The first assumption is the linear dependence of 
Vhxc on the occupation numbers 

VSxcK]=P^VSxcK■,P^^^]■ (6) 

Above procedure is exact for the Hartree part of the 
potential, but it is approximate for the much smaller 
exchange-correlation part. In this point the orbital 
exchange-correlation potential has to be calculated with 
fully spin polarized orbital density. 

2) The second simplification assumes employing the 
spherically averaged radial local orbital density nf (r) to 
compute the local orbital potential V§xc- 

Therefore, the angular part characterized by quantum 
number m; of pseudo-orbitals is used only to calculate pf 
and Cf as follows 

7/,,„„/(r) = i Pl^a VSxcinliir); 1] (r), 

C/V,,/ = I Plraa I dr VSxcKiAr); 1] (0?,™,^(r))^ 

where the indices to/ denote the angular momentum 
quantum number of the shell I (s, p, d, or /) of the atom 
of type /. 

The total energy within the non-variational spin po- 
larized pSIC scheme is constructed to resemble the DFT 
one and reads 

Esic[n, m]=^ /^k £nk + Eton - 

^ / drn^r)VSxcMr),m{r)]+EHxc[n{r),m{r)]- 

E EhXcK] + E (Ck|l>S/clCk>, (7) 

i.(T nk,o" 
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where n(r) and m(r) are the total and the spin polariza- 
tion density, respectively. 

The exchange-correlation part of the total energy cor- 
rection is a small number defined as 

EhxcK] ^ Jdr <(r) Qf^K^] + sxcKW]) , 

where Sxc is the local exchange-correlation energy den- 
sity. 

The last term in the formula ([7]) is the band correction, 
and shifts the total energy very strongly, restoring its 
proper curvature with respect to a change of the lattice 
constant (see Fig. 7 in Ref. [11]). 

In the scheme presented here, we implement equations 
for the ultrasoft pseudopotentialSfSl since they allow for 
substantial reduction of the energy cutoff for systems con- 
sisting of transition metals and rare earth atoms. How- 
ever, the USPP are not norm- conserving and need some 
additional terms to be included in the ordinary DFT and 
the pSIC methods. These terms contain the augmented 
charges Qaa' and projectors f3a. The overlap matrix for 
an orthonormality condition is 



S 



(8) 



where 



Qaa' ^ J '^^ Qaa'{r) 

(f)-^^ and (j)^^ are the all-electron and the pseudo-atomic 
functions, and a = [n,Z,m,/] sets all quantum numbers 
for the atom /. 

The pseudopotential splits into the local part Vloc (r) 
and the non-local part D"^^,, which consists of the 
Kleinman-Bylander term D'^^, and the augmentation 
term as follows 



Kc' + / dv (VLocir) + V^xcir)) Qao'(r). 



With the above definitions, the pSIC orbital density is 
<{v)=p1 ( |(/..(r)p + Q„„,(r) (/3„,|C> ), 

OLO.' 

and the pSIC-USPP occupation numbers are 

pr = E/"k(Cki</'.)(0.iCk) X 

nk 
nk 



The pSIC potential within the USPP scheme contains 
an additional term which reads 



i OiQ.' 



Thus, the Kohn-Sham equation with the USPP is 



The total energy terms of the pSIC origin are 



-E^^^c:K]+E/"k (Ckl iVs^ic + Vffs) ICk)- (9) 

2,cr nk 

In addition, the pSIC equations in the covariant form 
contain the off-diagonal occupation numbers 



Pi .nil inip/ 



V — 
^SIC — 



E/"k (V'^kl S <l>I.,mi.l){<l>I,m',U S* IV'J 
nk 

,1711 ,7n', ,1 \ '-^T^pM 



.k/j 



E 



^1/2 ^1/2 



I .mi ,mj 

7/,mi,; = Vnxcini.iir); 1] 0/,m,,/(r) 

Ci,mi,i = y dr 0/,„,,/(r) VHxc[ni,iir);l](t)i,mi,i{r). 

Another approximation for the augmentation part of 
the pSIC potential is made, assuming that the chosen 
pseudo-orbitals form a complete basis set. Thus, the 
beta projetors act on the atomic radial functions and en- 
able simple calculation of the radial integrals. Later, the 
Kohn-Sham states are projected onto the pseudo-orbitals 
in the plane-wave representation, as it is a case in the 
occupation numbers. The corresponding definitions are 
following 

VUS= E \Sh.miA)x 



9 Pl,mi,m[,l ^I,7n'mi'\l Wl..mi'\l ^ h 



and 



' 1 .7711 -.'^[i^ 



E {<I^I,mi.l\Pa) X 

OlOl' 
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In the above form, the V^g potential is computation- 
ally as simple as the occupation numbers, because the 
quantities e""^ depend only on the pseudopotential pa- 
rameters and can be calculated ones. 



III. FORCES IN THE PSIC SCHEME 

In this section, we give complete equations for forces 
in the pSIC scheme with ultrasoft pseudopotentials. 

According to the Hellman-Feynman theorem, the 
forces contain only the derivatives of the potentials and 
not the Bloch functions. The a index denotes one of 
the cartesian directions x,y, z from now on, and the a 
component of the displacement of atom / is denoted as 

Taj- 

Thus, following the equation ([SJ, we get an expression 
for the pSIC contribution to forces 



pSIC 



dEsic ^ V- ^^gxcKm,,;] 



dTa, 



/nk 

n,k 



OTaJ 



OTaJ 



The explicit derivatives are: 

dTaJ 

r dxf 

nil (T "/,m;.m;,/ r / \ i 

2 J drpj„^^^^,i ^^^^ X 
^Vff[n/,;(r); 1] + exc[ni,i{r); 1] 



(10) 



and 



1 ^-1/2 ^-1/2 

9 ^I,mi,l '-'I,m'l 



dVsic 

dTaJ 2 

ll,m,,l) 7, (7/ 



9r< 



aj 



1 ,7ni ,1 



dTa, 



Pi 



mi ,mj 



,/ + C.C. 



(11) 



and the ultrasoft part 
dVus 



dTa, 



_1 aug 

9 / y I,m'm,i", 



,771; " ,1,(7 



dpi I , 

\S(l)I,mi,l) ^ '' (<S'*0/,m,",/| -I- 



d{S4>I,mi,l) 



dTa, I 



dTa 



P'l,mi,m[,l ('S'*0/,,77,",/| + C.C. 



(12) 



The 



derivatives 



dpl-m^m'^il^^o^j and 

\d{S(f)i^m,j)/dTa,i) are defined in Ref. [21] by eqs. (13- 
19), and we give them explicitely in the appendix. The 
derivative \dji,mi,i/dTa,i) is obtained in the same way as 
the derivative \d{S4>i ,mi.i) / dTa,i) , because the potential 
VHXcV''i,i{f)\ moves together with the atomic functions. 

For the derivative of the overlap operator S, the follow- 
ing approximation is made. It is assumed that contribu- 
tions of the beta functions centred at the atoms different 
than the moved atom are neglected. It turns out that this 
approximation does not corrupt the accuracy, and it will 
be shown in the test cases later on. This simplification 
is necessary, because in the pSIC scheme the projectors 
used in the definition of the occupation numbers have 
to be orthogonalized, which in turn sets a difficulty in 
calculation of the derivatives. 

Above definitions are valid for the non- variational pSIC 
approach. First approximate equations for forces have 
been given by Filippetti and Fiorentini^ however their 
formalae neglected terms with the derivatives of the oc- 
cupation numbers. Recent work by Filippetti et al}^ for 
the variational pSIC scheme contains similar expressions 
for forces. We have added the derivatives of occupation 
numbers in a way akin to the equations for forces in the 
DFT+U schemCfSi These terms are rather small, and we 
show their effect discussing the Ce02 case in the next 
section. 



IV. TESTS FOR FORCES AND RELAXATIONS 



Wurzite ZnO and CeOo in the F„ 



structure 



As a first test case, we employ introduced scheme for 
forces to the wurzite ZnO. We use the ultrasoft pseu- 
dopotential, the LDA exchange-correlation functional in 
the parametrization of Perdew-Zunger, the kinetic energy 
cutoff of 35 Ry, and the uniform Monkhorst-Pack (6,6,6) 
k-mesh in these calculations. 

The results for the total energy and the force acting 
on the displaced atom Zn(l) in the wurtzite unit cell are 
presented in Figure [TJ The Zn atom is displaced only in 
the z-direction and the magnitude of the displacement 
is given as a function of Wyckoff position in units of the 
lattice constant. 

First, we discuss the role of the approximation sim- 
plifying the orthogonalization of local atomic projectors 
on the total energy vs. atomic displacement curves and 
forces for both LDA-I-U and pSIC methods. As we have 
mentioned in the section IIIII in this approximation the 
non-local contributions of beta functions to the deriva- 
tives are neglected, and only the diagonal terms in the 
beta functions are considered when the derivative with 
respect to the atomic position is calculated. The LDA+U 
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calculations (with U=5 eV) with non-orthogonalized pro- 
jectors, called "atomic", are performed without any ap- 
proximation. Simultaneously, calculations of the approx- 
imate forces obtained with the orthogonalized projectors, 
called " ortho- atomic" , are compared to results from the 
exact formulae. Panels a) and b) of Figure [1] show a 
perfect agreement between the results for the two sets of 
projectors applied for the c?-shell, ensuring us that the ap- 
plied approximation for the derivatives in forces is rather 
good. 

In panels c) and d) of Figure [TJ the pSIC results are 
presented for the same atomic displacements which have 
been described above for the LDA-I-U method. As one 
can see, the force vanishes exactly at the geometry that 
coincides with the atomic position for which the total 
energy gets the minimum. It clearly demonstrates the 
correctness of the equations for forces derived for the 
pSIC method in this paper. 



Exp. 



LDA LDA+U pSIC 



lattice constant a 6.16 6.04 
lattice constant c 9.84 9.68 

starting distorted geometry 

Zn(l)-0(1) - 3.382 

Zn(2)-0(2) - 3.986 

relaxed parameters 

u parameter 0.382 0.381 

Zn-O bond 3.759 3.684 



6.05 
9.70 



3.388 
3.993 



0.381 
3.691 



6.09 

9.76 



3.410 
4.019 



0.379 
3.697 



TABLE I. The geometry parameters of initially distorted 
and fully relaxed wurzite ZnO, calculated with the LDA, the 
LDA+U and the pSIC methods; obtained with the BEGS al- 
gorithm which contains forces. Lattice constants and bond 
lengths are in a.u. The experimental values are from Ref. 
[25]. 



Next, a relaxation of the displaced atomic positions 
within the wurzite ZnO cell has been performed within 
the Newton-Raphson optimization scheme based on the 
Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm^i 
for the estimate of the inverse Hessian matrix. The cri- 
teria for the geometry optimization have been set as: 
the energy difference between subsequent BFGS-steps < 
10-4 Ry, and the force < IQ-^ Ry/a.u. 

The starting non-equilibrium geometry has been ob- 
tained by application of the same distortion for all cal- 
culations: the LDA, the LDA-I-U, and the pSIC. The 
wurzite structure has been perturbed in a such way that 
two "u" parameters (that determine the Zn-O distances 
along the z-axis and are defined as the bond length along 
the hexagonal symmetry axis devised by the lattice con- 
stant c) for the Wyckoff positions have been chosen for 
the Zn-O distances along the z-axis: to be equal to 
wi=0.349 and M2=0.412, which corresponds to consid- 
erably shorter and longer bond lengths, respectively (for 
the perfect wurzite structure ui=U2). The lattice con- 
stant has been optimized for each method prior to the 
relaxation. The identical initial distorted geometry has 
been used to find the equilibrium geometry within the 
LDA, LDA-I-U, and pSIC approaches. 

The results are displaced in Table U that collects the 
optimized lattice constants, the Zn-O bond lengths along 
the z-axis for the distorted and relaxed structures, and 
the optimized u parameters. For all three methods, the 
optimized u parameters (i.e., ui and U2) are identical. 
This correct result strongly corroborates the correctness 
of the derived equations for forces in the pSIC scheme. 
Note also that the lattice constant c obtained in pSIC 
method agrees better with the experimental value (9.83 
a.u.) than the lattice constants obtained in the LDA 
and the LDA-I-U schemes. 

As the next test case for the equations for forces, 
we consider a rare earth compound Ce02 in the FmSm 
structure. Here, we have chosen for the calculations 



the USPP, the Perdew-Zunger LDA functional, the ki- 
netic energy cutoff of 35 Ry, and the uniform (8,8,8) 
Monkhorst-Pack k-mesh. 

In Figure [2l the total energy and the force acting 
on the displaced atom Ce(l) are shown. The atom is 
displaced only along the [1,1,0] crystal- axis and the 
magnitude of the displacement is given as a function 
of Wyckoff position in units of the lattice constant, 
which has been fixed for this test at the experimental 
value of 5.41 a.u. The total energy minimum and the 
zero force occur exactly at the equilibrium geometry 
of the non-distorted structure. This is a next proof 
for the derived force formulae, which work also for 
the /-electron compound. For a comparison, we have 
also presented approximate forces which have been 
obtained neglecting the derivatives with respect to the 
atomic position of the occupation numbers. Such terms 
enter equations (llOH12p . and they were omitted in 
Refs. [12,18]. However, similar terms are present in the 
DFT-I-U forces. In these two cases, the forces slightly 
differ, but both approaches give zero force at the same 
geometry. 



B. Distorted perovskites: LaTiOs, YTiOs, LaMnOs 

Strongly correlated perovskites LaTiOa, YTiOa, and 
LaMnOa exhibit Jahn- Teller distorsions and crystallize 
in the Pnma structure. They have been widely studied 
within the DFT-f U method, just to mention as an exam- 
ple the work by Okatov et alM (for LaTiOa and YTiOs), 
and by Trimarchi and Binggeli^^ (for LaMnOs). 

Nevertheless, the self-interaction correction applied to 
the oxygen atom in these compounds may cause some 
changes in the predicted geometry in comparison to the 
DFT-hU resuhs. 

At low temperatures, LaTiOa has a G-type antiferro- 
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Atom Class 



Coordinated 



RE,Oi 4c 
TM 4a 
O2 8d 



{u,l/A,v), (u+l/2,3/4,u+l/2) 
{U,Z/A,v), (it+l/2,l/4,t;+l/2) 

(0,0,0), (1/2,0,1/2) 
(0,1/2,0), (1/2,1/2,1/2) 

±{x,y,z) 
±(a;,y,2) + (l/2,0,l/2) 
±(a;,3/,2)+(0,l/2,0) 
±(x,y,2)+(l/2,l/2,l/2) 



TABLE II. The Wyckofl positions for each ionic specie in the 
Pnma structure of RETMO3 (RE=La,Y and TM=Ti,Mn). 



magnetic structure and YTiOa is a ferromagnet, while 
a colossal magnetoresistance material LaMnOs is an A- 
type antiferromagnet. It is known that relations between 
the cell-axes determine the magnetic order in distorted 
perovskites. However, the calculation of stress tensor 
is not implemented yet in the pSIC approach. There- 
fore, we focus on the FM-ordered structures keeping the 
cell parameters fixed at the room-temperature crystallo- 
graphic data. The details of the Pnma crystal structure 
are given in the Table |lTl For such cell, we optimized the 
geometry employing various DFT schemes, namely the 
GGA, the GGA-HU, and the pSIC. 

For all schemes, we have chosen the Perdew-Burke- 
Ernzerhof functional and employed the ultrasoft pseu- 
dopotentials. In the case of the GGA-I-U, the Hubbard-U 
parameter for Ti and Mn was set to 3.0 eV. In the pSIC 
calculations, the self-interaction correction has been ap- 
plied to the outermost d-shell of rare-earth (RE) and 
transition-metal (TM) elements and also to the 2s- and 
2p-shells of the oxygen. The results of calculation within 
the GGA, the GGA+U, and the pSIC methods are col- 
lected in Table IIIIl which presents crystallographic pa- 
rameters obtained from the BFGS optimization and com- 
pares them to the experimental data. 

As one can see, the distortions calculated with the 
pSIC method are usually larger than obtained from the 
GGA and the GGA+U methods. Most of structural pa- 
rameters calculated within the pSIC method are closer 
to the GGA-I-U numbers than to the GGA ones. Nev- 
ertheless, inclusion of the self-interaction correction to 
the 2s and 2p shells of the oxygen leads to a substantial 
difference, and brings the pSIC results usually closer to 
the experimental values. Some discrepancies still exist, 
especially for small parameters, and their reasons may 
lay on the accuracy of either the theoretical methods or 
experimental techniques. On the theoretical side, for ex- 
ample, the reported calculations involve the pseudopo- 
tentials and it is an open question how obtained results 
would differ from the results of all-electron approach. 

Table IIIIl gives also the mean error of the calculated 
parameters P^^ with respect to the experimental values 



Parameters Exp. 


GGA 


GGA+U 


pSIC 


LaTiOS 












Exp.^ a= 


10.6647 a.u., b= 


=14.9300 


a.u. 


c=10.5607 


a.u. 


RE u 


0.4916 


0.4635 




0.4685 


0.4734 


RE V 


0.0457 


-0.0014 




-0.0019 


0.0107 


Ol u 


0.0799 


0.0163 




0.0330 


0.0609 


Oi V 


0.0087 


-0 0464 




-0.0862 


0.0877 


O2 X 


0.2096 


0.2022 




0.1938 


0.1997 


O2 y 


0.0417 


0.0259 




0.0400 


0.0374 


O2 z 


0.2941 


0.2920 




0.3124 


0.3317 


YTi03 












Exp.i^ a= 


10.0375 a.u., b= 


=14.3827 


a.u. 


c=10.7318 


a.u. 


RE u 


0.4793 


0.4486 




0.4317 


0.4633 


RE V 


0.0729 


-0.0076 




0.0131 


0.0109 


Ol u 


0.1211 


0.0268 




0.0266 


0.0558 


Ol V 


0.0042 


-0 0980 




-0.1227 


0.0921 


O2 X 


0.1910 


0.1852 




0.1864 


0.1791 


O2 y 


0.0580 


0.0470 




0.0642 


0.0358 


O2 z 


0.3100 


0.3114 




0.3062 


0.3449 


LaMnOS 












Exp.^^ a= 


10.8508 a.u., b= 


= 14.4904 


a.u. 


c=10.4540 


a.u. 


RE u 


0.5490 


0.5525 




0.5536 


0.5524 


RE V 


0.0100 


0.0097 




0.0101 


0.0093 


Ol u 


-0.0140 


-0.0211 




-0.0232 


0.0255 


Ol V 


-0.0700 


-0.0834 




-0.0910 


0.0791 


O2 X 


0.3090 


0.2990 




0.3068 


0.3192 


O2 y 


0.0390 


0.0434 




0.0458 


0.0436 


O2 z 


0.2240 


0.2144 




0.2180 


0.2274 


Mean error (A) 










RE u 




0.13 




0.16 


0.08 


RE V 




2.24 




1.65 


2.37 


Ol u 




2.08 




2.02 


1.60 


Ol V 




30.86 




41.42 


34.14 


O2 X 




0.10 




0.11 


0.14 


O2 y 




0.68 




0.32 


0.60 


O2 z 




0.05 




0.10 


0.25 



TABLE III. Experimental and theoretical parameters of 
Pnma structure (in crystal coordinates) for LaTiOs, YTiOs 
and LaMnOs, obtained from the BFGS relaxation within the 
GGA, the GGA+U, and the pSIC methods. Last block in the 
table gives the mean error of parameters obtained with each 
method in comparison to the experimental values (defined by 
eq. (HSI). 



pExp.^ it is defined as 

(A>= E 



P 



M 



pExp 



pExp 



(13) 



where the summation runs over all calculated structures: 
LaTiOs, YTiOs, LaMnOa. One general observation is 
clear: the smaller is the parameter, the larger is the dis- 
crepancy between the calculated and the experimental 
values. Generally the distorsions from the ideal per- 
ovskite structure are larger in the calculations than in 
the experiment. This might be due to the fact that, in 
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111 




220 






TM-TM 


TM-Si 


TM-TM 


TM-Si 


ideal Si geom. 


4.4686 


4.4686 


7.2983 


4.4686 


GGA, TM= 




4.8132 


4.5202 


7.1868 


4.4871 


pSIC, TM= 


=Mn 


4.9495 


4.5429 


7.2735 


4.5552 


GGA, TM= 


=Re 


4.1837 


4.5532 


5.7524 


4.4768 


pSIC, TM= 


=Re 


4.1713 


4.5181 


6.5140 


4.4108 



TABLE IV. The distances TM-TM and TM-Si (in Bohr) in 
Si:Mn and Si:Re for two configurations of impurities: 111 and 
220 obtained after the BFGS minimization from the GGA 
and the pSIC approaches. 



the experiment, the signal is averaged over the sample, 
which is never clean and so ideally periodic like in the 
calculations. 

Concerning the FM-order, all theoretical methods 
give the magnetic moments of the Ti atom equal to 1.0 
/is in both LaTiOs and YTiOs, whereas the calculated 
magnetic moment at Mn in LaMnOa is 4.0 fj,B- Discus- 
sion of magnetic structure issues runs beyond the scope 
of this work, however, we would like to mention that 
the results obtained in this paper agree with numbers 
calculated within the GGA and the GGA-f-U schemes 
and reported earlier by other authorsi^ 



Diluted magnetic semiconductors: 
Si:Mn and Si:Re 



TM-Si-TM (for 220); obtained from the BFGS minimiza- 
tion performed in the GGA and the pSIC schemes, and 
compared to the ideal geometry of the silicon crystal. 

In the case of the close distance pairs (111), the Mn 
ions repel themselves, while the Re ions attract each 
other in comparison to distances in the ideal silicon crys- 
tal. This effect is considerably stronger in the pSIC than 
in the GGA method. 

For the 220 pairs, the TM ions get closer in the both 
cases of Mn-Mn and Re-Re pairs, the effect being espe- 
cially pronounced for Re ions. In contrast to the 111 case, 
this attraction of TM pairs effect is much weaker in the 
pSIC than in the GGA approach. The TM-Si distances 
usually become slightly longer than the ideal Si- Si bond, 
except for the Re-Si-Re bridge in the pSIC approach. 
This effect is important for the magnetic properties of 
silicon doped with Re and will be published elsewhere. 
Here, we only comment on the fact that, the rhenium ions 
in silicon have smaller magnetic moment (1 /i^) than the 
Mn ions (3 /is), and therefore, rhenium employs more 
valence electrons for a hybridization with atoms of the 
host and with another close Re ion. Due to a larger lo- 
calization of the d-shell electrons in Re within the pSIC 
approach, these states contribute much weaker to a hy- 
bridization between Re- Re, and this bond is much longer 
than in the GGA method. A very interesting difference 
between Si:Mn and Si:Re is in the DOS: the states, which 
are closer to the Fermi level, originate from the closest 
neighbours of the impurity in the case of Mn, and from 
the second close neighbours in the case of Re. This fact 
gives one of the reasons why the 220 pair of Re in Si 
relaxes stronger than the 111 pair. 



As the third example, we have chosen two prototypes 
of the DMS systems. We consider the silicon crystal 
doped (i) with two Mn, and (ii) two Re impurities per 
cell. Detailed investigations of structural and magnetic 
properties of these DMS's will be given elsewhere. 
Here, we only present an effect of the pSIC scheme on 
the geometry around the transition-metal ions (TM) 
by comparing the atomic positions obtained from the 
pSIC and the standard GGA method. We consider two 
geometries of the TM pairs substituted into Si sites 
within the cubic unit cell with 64 atoms (with the silicon 
lattice constant resulting from the GGA calculations 
and equal to 10.32 a.u.). We consider (i) two TM atoms 
being the nearest neighbors (hereafter indicated as 111, 
since they take the sites (000) and a/4(lll) in the silicon 
crystal, where a is the silicon lattice constant) and (ii) 
two TM atoms in the next nearest neighbours sites, they 
are bridged by the Si atom (hereafter indicated as 220, 
since they occupy the sites (000) and a/4(220)). 

Table lTVl presents the distances between: (i) transition- 
metal ions (TM-TM), and (ii) the transition metal and 
the silicon atom adjacent to the one of the TM-ions (for 
111), and (iii) the TM-ion and the Si atom at the bridge 



V. SUMMARY 

We have derived the expressions for forces within the 
non- variational pSIC approach with ultrasoft pseudopo- 
tentials used to account for electron and ion interac- 
tions and implemented the scheme into the quantum 
ESPRESSO plane-wave code. First, we have performed 
benchmark calculations to check the internal consistency 
of the scheme for the wurzite ZnO and rare-earth /- 
electron compound Ce02 in the FmSm structure. In both 
cases, the forces within the pSIC scheme vanish for the 
geometry corresponding to the minimum of the total en- 
ergy. Also optimization procedure within the code works 
perfectly bringing the initially distorted crystallographic 
structures of ZnO and Ce02 into the correct equilibrium 
geometry efficiently. 

Further, we have performed calculation within the 
pSIC approach to determine the geometry of distorted 
perovskites LaTiOa, YTiOa, and colossal magnetoresis- 
tance compound LaMnOa in the Pmna structure, and also 
of silicon doped with pairs of Mn and Re ions. These sys- 
tems have been chosen, since there are indications that 
the spurious self-interaction and resulting more diffused 
electronic states can lead to certain systematic errors. 
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Indeed, in the cases studied here, the pSIC results for ge- 
ometry parameters are usually closer to the experimen- 
tal ones than the parameters obtained from the stan- 
dard approximations of the DFT and the DFT-I-U meth- 
ods. This strongly suggests that the larger localization of 
the electronic states is better accounted for in the pSIC 
scheme, which could provide also more reliable predic- 
tions in many systems. Also in the case of Mn and Re 
pairs in silicon, the geometries of the systems obtained 
within the pSIC and the GGA differ considerably. Ef- 
fect of the pSIC relaxations is usually weaker than the 
GGA ones, which is a consequence of weaker sp-d hy- 
bridization. An exception is the Re-Si-Re configuration 
for which the Rel(5d)-Re2(5d) interactions are strong 
and the pSIC relaxations are larger than those obtained 
from the GGA method. 

Having functioning scheme to calculate forces within 
the pSIC method, the further studies are under way to 
determine the areas of relevant applications and deeper 
investigate the reliability of the method. 
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Appendix A: Derivatives of the orbital occupation 
numbers with respect to the ionic displacement 

Partial derivatives of the occupation numbers. 
Pi mi m' J' ■^ith respect to the atomic displacements, Tq /, 
are given in Ref. [21] by eqs. (13-19). Nevertheless, for 
the completeness, we include these derivations here. 

We start from the occupation numbers in the norm- 
conserving pseudopotential scheme. 



nil ^nii ,/ 



dTaJ 
fnk 



dTaJ 



((V'nkl0/,m,,i>)(</'/,m,,i|^„k/ 



d 



,mi,l) ^ {4'I,mi,l 
OTa.I 



The derivative of {'ipnk\'t'i,mi.i} reduces to the derivative 
of 4>i,mi,i, since due to Hellman-Feynman theorem ■f/'J^^. 
does not change with the displacement. 



The atomic orbitals 4>i,mi,i are represented in the 
plane- wave basis at each vector k from the IBZ, in or- 
der to project them onto the Bloch functions. Then, the 
projection is symmetrized, to take care of the summation 
over all points from the BZ. The atomic orbital at point 
k is expressed: 

0/,m,,i.k(r) = ^^e"*'''^0/,„,,;(r-R-Ti) = 



R 



-zk-r 



1 ^^-.k(r-R)^^^^_^^(^_j^„^^)^ 



R 

N is the number of the direct lattice vectors R. The func- 
tion (/'/,m,./(r — R. — '''i) is periodic with the lattice and 
its Fourier expansion in the reciprocal lattice vectors G 
is defined as: 



<t>hmi,i{v) 



—y 



e-'(*'+°)-'-cz,„„,(k + G), 



where V is the volume of the system (V^Nfi and 51 is 
the cell volume). The Fourier components _;(k+ G) 
read: 

c/,„,,(k + G) = y dr e'(k+G)T ^^^^^ ^(r) = 



1 



„»(k+G)-ri 



R 



„i:(k+G)-r 



'/'/,m,j(r)- 



The derivative of the atomic function (i>i^mi.i with respect 
to the displacement of the same atom I in the direction 
a is thus 



m; ,1 



where (k + G)^ is the vector component along the polar- 
ization a. 

The derivatives of the occupation numbers in the 
norm-conserving pseudopotential scheme are nonvanish- 
ing only for the displacement of the same atom at which 
the occupations are considered. 

In the ultrasoft-pseudopotential scheme, the deriva- 
tives \d{S(l)i,mi,i)/dTaj) have to be computed. Accord- 
ing to eq. ([U, the above derivative contains derivatives 
of the Pa functions (here the index a = [n,l,m, I]). 
These functions are the ultrasoft pseudopotential func- 
tions, which can be expressed also in the plane-wave rep- 
resentation. The overlap given by eq. ^ is nonlocal in 
/3a- Therefore, we made the approximation mentioned 
in section III, and we neglected contributions from the 
derivatives of the /3a' functions centred at atoms I' dif- 
ferent than the moved atom I. 
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FIG. 1. Total energy of the wurzite ZnO and the force acting 
on the displaced atom Zn(l) along z-axis. Panels a) and b) 
are for the LDA+U (U=5 eV) with "atomic" and "orthogo- 
nalized atomic" projectors for USPP, panels c) and d) are for 
the LDA+pSlC with " orthogonalized atomic" projectors for 
USPP and NCPP. Dotted lines are just guides for the eye. 
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FIG. 2. Total energy of Ce02 in the Fm3m structure and 
the forces (according to our pSIC equations - squares, and 
approximated witliout terms containing the derivatives of the 
occupation numbers - triangles) acting on the atom Ce(l) 
displaced along the [1,1,0] crystal-direction; calculated with 
the USPP and the orthonormalized projectors. 



